** Set filepaths to folders to store intermediate files, output and .do files:
global datafiles "N:\SpecialProjects\Inequality_Egypt\Peru_2019\data"
global dofiles "N:\SpecialProjects\Inequality_Egypt\Peru_2019\do_files"

* The following .do file runs the estimation code for DHSV and populates a matrix with the estimates (before producing bootstrapped standard errors).  These are the estimates found in Table 3 (and B.3 and B.4).
do "$dofiles\DHIS_IV_pre_boot_EDCC_replication_file.do"


use "$datafiles/peru_dhsIV_analysis.dta", clear


**** The bootstrap routine starts here ****


/*** Set up a program called myboot to run the estimation for the bootstrap ***/
capture program drop myboot
	program define myboot, rclass
	preserve
	bsample
	  
	********************* Now execute the estimation procedure *********************************  
	
	** Policy and analysis years for main specification (as seen in Table 3)
	local aby "1990"
	local aey "1997"
	local pby "1996"
	local pey "1997"
	local exc "& year~=1995"

	/*
	** alternative specifications (Table 6):

	** extend policy years from 1995 to 1998:
	local aby "1990"
	local aey "1998"
	local pby "1995"
	local pey "1998"
	local exc " "

	** extend policy years from 1995 to 1999:
	local aby "1990"
	local aey "1999"
	local pby "1995"
	local pey "1999"
	local exc " "

	*/

		gen bid=_n

		tempname bootsample
		save `bootsample', replace

	************* Logit (as described in Section 6 and shown in Table 1) ******************************
	keep bid age_* ster_* kids_* mar_* will_end*  boys_*  born_* m_age_fb yr_st ed_cat weight1 rur_coast urb_mount rur_mount urb_jungle rur_jungle urb_coast 

	*Create the pseudo-panel
	reshape long age_ ster_ kids_  mar_ will_end_  boys_  born_  , i(bid) j(year)

	* dummy for being sterilized in a given year
	gen sterilized=yr_st==year


	* Drop observations after sterilization (no longer at risk of sterilization)
	bys bid: gen flag=cond(sterilized==0, sum(sterilized),0)
	drop if flag==1


	* It was "illegal" to be sterilized if you were not married and did not have kids - and almost noone was, so we make this 
	* sample restriction (otherwise the probit does not converge due to perfect prediction of failure.)
	gen eligable=(mar_ | will_end_) & kids_>0


	gen kids_count=kids_
	replace kids_count=. if kids_==0
	replace kids_count=10 if kids_>10
	tab kids_count, gen(kids)

	gen year_lin=year-1992

	gen young=age_<26
	gen mid1_age=age_>=26 & age_<=29
	gen mid2_age=age_>29 & age_<=36
	gen older=age_>36

	gen ed_prim=ed_cat==2 | ed_cat==1
	gen ed_second=ed_cat==3
	gen ed_high=ed_cat==4

	gen no_birth=born_==0

	gen no_boys=boys_==0


	local int_vars " rur_mount rur_jungle ed_prim young mid1_age no_birth "

	cap drop policy
	gen policy=0
	replace policy=1 if year>=`pby' & year<=`pey'


	cap drop policy_X_*
	foreach x of local int_vars{
	qui gen policy_X_`x'=`x'*policy
	}

	forval i=1/9{
	gen policy_X_kids`i'=policy*kids`i'
	}


	gen policy_save=.
	foreach x of local int_vars{
	qui gen policy_X_`x'_save=.
	}

	forval i=1/9{
	gen policy_X_kids`i'_save=.
	}


	logit sterilized  kids1-kids9 m_age_fb rur_coast urb_mount rur_mount urb_jungle rur_jungle ed_prim ed_second young mid1_age older  no_boys no_birth  ///
		  policy_X_rur_mount policy_X_rur_jungle policy_X_ed_prim policy_X_young policy_X_mid1_age  policy_X_no_birth  year_lin [pw=weight1] if year>=`aby' & year<=`aey' `exc' & eligable, vce(cluster bid)


	qui cap drop pr2_

	qui predict pr2_ if year>=`pby' & year<=`pey' `exc' & eligable & e(sample)


	qui replace policy_save=policy
	foreach x of local int_vars{
	qui replace policy_X_`x'_save=policy_X_`x'
	}

	forval i=1/9{
	qui replace policy_X_kids`i'_save=policy_X_kids`i'
	}


	qui replace policy=0
	foreach x of local int_vars{
	qui replace policy_X_`x'=0
	}

	forval i=1/9{
	qui replace policy_X_kids`i'=0
	}

	qui cap drop pr1_


	qui predict pr1_ if year>=`pby' & year<=`pey' `exc'  & eligable & e(sample)


	qui replace policy=policy_save
	foreach x of local int_vars{
	qui replace policy_X_`x'=policy_X_`x'_save
	}

	forval i=1/9{
	qui replace policy_X_kids`i'=policy_X_kids`i'_save
	}

	  
	save "$datafiles/peru_dhsIV_qp.dta", replace
	
** All de jure women who were eiligable in 1996 or 1997 have a pr1 pr2 for that year 
** Being eligible means having ever been married and having at least one child	

	keep pr1_* pr2_* bid year 
	
	* make this wide again
	reshape wide pr1_ pr2_  , i(bid) j(year)
	sort bid
	save  "$datafiles/pscoresIV.dta", replace

	** Merge the p1 p2 variables onto the cross-sectional data set.
	use `bootsample', clear
	cap drop _merge
	sort bid
	merge bid using "$datafiles/pscoresIV.dta"

	cap drop S
	gen S = yr_st>=`pby' & yr_st<=`pey' 


	gen status=.
	replace status=1 if yr_st<`pby' & yr_st~=.               /* sterilized before the policy -- these women are all excluded from the outcome results below */


	egen _p2=rowmean(pr2_*)
	egen _p1=rowmean(pr1_*)

	gen delta_p=_p2-_p1
	gen dp_p2=delta_p/_p2

	*** Standard Reweighting 

	gen sdrweight=weight1 if S==1 & status~=1
	replace sdrweight=weight1*(_p2/(1-_p2)) if S==0 & status~=1

	bys S: egen sdsrweight_t=total(sdrweight)
	gen sdsrweight=sdrweight/sdsrweight_t


	*** Our reweighting S2=1, D=1
	gen rweight=weight1*dp_p2 if S==1 & status~=1                 
	replace rweight=weight1*(delta_p/(1-_p2)) if S==0 & status~=1 

	bys S: egen srweight_t=total(rweight)
	gen srweight=rweight/srweight_t

	*** Our reweighting S2=1, D=0

	gen orweight=weight1*(_p1/_p2) if S==1 & status~=1           
	replace orweight=weight1*((_p1)/(1-_p2)) if S==0 & status~=1 

	bys S: egen osrweight_t=total(orweight)
	gen osrweight=orweight/osrweight_t



	tempname bootsample_ww
	save `bootsample_ww', replace
	
	
	******** estimating the number affected *********************
	gen w3=weight1*246.9

	total w3 if S==1
	matrix total=e(b)
	return scalar women_all=total[1,1]

	gen estim=w3*dp_p2
	total estim if S==1
	return scalar women_campaign=_b[estim]


	**************************** Number of Kids *****************************************************************
	* Modified Reweighting (S2=1, C=1) 
	reg num_kid S  [pw=srweight] if status~=1 
	return scalar kids_s2d1=_b[S]

	* Modified Reweighting (S2=1, C=0) - only positive weigths are included by default since pw drops negative weights
	reg num_kid S [pw=osrweight] if status~=1 
	return scalar kids_s2d0=_b[S]



	**************************** Mother's Labor-Force Participation ((working for pay) *******************************


	* Modified Reweighting (S2=1, C=1) 
	reg working_p S  [pw=srweight] if status~=1 
	return scalar wp_s2d1=_b[S]


	* Modified Reweighting (S2=1, C=0) - only positive weigths are included by default since pw drops negative weights
	reg working_p S [pw=osrweight] if status~=1 
	return scalar wp_s2d0=_b[S]


	**************************** Mother's Health *****************************************************************
				 ** weight_height_m

	* Modified Reweighting (S2=1, C=1) 
	reg weight_height_m S  [pw=srweight] if status~=1 & pregnant==0
	return scalar WH_M_s2d1=_b[S]

	* Modified Reweighting (S2=1, C=0) - only positive weigths are included by default since pw drops negative weights
	reg weight_height_m S [pw=osrweight] if status~=1  & pregnant==0
	return scalar WH_M_s2d0=_b[S]

				 ** height_age_m

	* Modified Reweighting (S2=1, C=1) 
	reg height_age_m S  [pw=srweight] if status~=1 & pregnant==0
	return scalar HA_M_s2d1=_b[S]

	* Modified Reweighting (S2=1, C=0) - only positive weigths are included by default since pw drops negative weights
	reg height_age_m S [pw=osrweight] if status~=1  & pregnant==0
	return scalar HA_M_s2d0=_b[S]

				   **Anemic

	* Modified Reweighting (S2=1, C=1) 
	reg anemic_m S  [pw=srweight] if status~=1 & pregnant==0
	return scalar anemic_M_s2d1=_b[S]

	* Modified Reweighting (S2=1, C=0) - only positive weigths are included by default since pw drops negative weights
	reg anemic_m S [pw=osrweight] if status~=1 & pregnant==0
	return scalar anemic_M_s2d0=_b[S]





	*********** Biometrics  W-f-H and H-f-A **********************************************************************
	** height for weight only recorded for kids 4 and under and we only want to look at kids born before the policy.  
	** So we can only do this for DHSIV and for these outcome the unit of analysis is the child (whose mother may or may
	** not have been sterilized) 

	
	drop kby_07-kby_14

	keep bid status srweight weight1 osrweight S  kby_*  kid_boy* kid_age* anemic_k* weight_height* height_age*

	reshape long  kby_  kid_boy kid_age anemic_k weight_height height_age, i(bid) j(person)

	keep if kid_age~=.

	*********************************** Weight for Height for kids *******************
	* Modified Reweighting (S2=1, C=1)
	reg weight_height S i.kid_age [pw=srweight]   if status~=1 & kby_<=1997 
	return scalar WHkid_s2d1=_b[S]


	* Modified Reweighting (S2=1, C=0) 
	reg weight_height S i.kid_age  [pw=osrweight]  if status~=1 & kby_<=1997 
	return scalar WHkid_s2d0=_b[S]


	******************Height for Age for kids**********************
	* Modified Reweighting (S2=1, C=1)
	reg height_age S i.kid_age [pw=srweight]   if status~=1 & kby_<=1997 
	return scalar HAkid_s2d1=_b[S]


	* Modified Reweighting (S2=1, C=0) 
	reg height_age S i.kid_age  [pw=osrweight]  if status~=1 & kby_<=1997 
	return scalar HAkid_s2d0=_b[S]



	******************Anemia for Kids**********************
	* Modified Reweighting (S2=1, C=1)
	reg anemic_k S i.kid_age [pw=srweight]   if status~=1 & kby_<=1997 
	return scalar Anemiakid_s2d1=_b[S]


	* Modified Reweighting (S2=1, C=0) 
	reg anemic_k S i.kid_age  [pw=osrweight]  if status~=1 & kby_<=1997 
	return scalar Anemiakid_s2d0=_b[S]



	******************Height for Age for girls**********************
	* Modified Reweighting (S2=1, C=1)
	reg height_age S i.kid_age [pw=srweight]   if status~=1 & kby_<=1997 & kid_boy==0 
	return scalar HAgirl_s2d1=_b[S]


	* Modified Reweighting (S2=1, C=0) 
	reg height_age S i.kid_age  [pw=osrweight]  if status~=1 & kby_<=1997 & kid_boy==0 
	return scalar HAgirl_s2d0=_b[S]



	******************Height for Age for boys*****************************************

	* Modified Reweighting (S2=1, C=1) -
	reg height_age S i.kid_age [pw=srweight] if status~=1 & kby_<=1997 & kid_boy==1
	return scalar HAboy_s2d1=_b[S]



	* Modified Reweighting (S2=1, C=0) 
	reg height_age S i.kid_age  [pw=osrweight] if status~=1 & kby_<=1997 & kid_boy==1 
	return scalar HAboy_s2d0=_b[S]



	************************************ Weight for Height for girls *******************
	* Modified Reweighting (S2=1, C=1)
	reg weight_height S i.kid_age [pw=srweight]   if status~=1 & kby_<=1997 & kid_boy==0 
	return scalar WHgirl_s2d1=_b[S]

	* Modified Reweighting (S2=1, C=0) 
	reg weight_height S i.kid_age  [pw=osrweight]  if status~=1 & kby_<=1997 & kid_boy==0 
	return scalar WHgirl_s2d0=_b[S]

	********************************** Weight for Height for boys *******************

	* Modified Reweighting (S2=1, C=1) 
	reg weight_height S i.kid_age [pw=srweight] if status~=1 & kby_<=1997 & kid_boy==1
	return scalar WHboy_s2d1=_b[S]

	* Modified Reweighting (S2=1, C=0) 
	reg weight_height S i.kid_age  [pw=osrweight] if status~=1 & kby_<=1997 & kid_boy==1 
	return scalar WHboy_s2d0=_b[S]


	***********EDUCATION OF CHILDREN

	**********Education of household children (under 15) but born prior to the program
	** need to use the household record to get these kids' education which leads to slight differences


	use `bootsample_ww', clear

	
	keep kid_age_ed* kbyed* kid_boy_ed* enrolled* sch_yr_* grade_at_age* bid status S srweight osrweight

	reshape long kid_age_ed kbyed kid_boy_ed enrolled sch_yr_ grade_at_age, i(bid) j(kid)

	rename kid_boy_ed kid_boy
	rename kid_age_ed kid_age
	rename kbyed kby




	*********** Education of own children**********************************************************************


	*********** Enrollment of own kids**********************************************************************

	* Modified Reweighting (S2=1, C=1) 
	reg enrolled S i.kid_age  [pw=srweight] if status~=1 & kby<=1997  & kid_age<=15 & kid_age>=6 
	return scalar ENkid_s2d1=_b[S]


	* Modified Reweighting (S2=1, C=0) 
	reg enrolled S i.kid_age [pw=osrweight] if status~=1 & kby<=1997  & kid_age<=15 & kid_age>=6 
	return scalar ENkid_s2d0=_b[S]

	gen ed_flag=e(sample)


	*********** grade for age of own children**********************************************************************
	
	* Modified Reweighting (S2=1, C=1) 
	qui reg grade_at_age S i.kid_age  [pw=srweight] if status~=1 & kby<=1997 & kid_age<=15 & kid_age>=6 & ed_flag==1
	return scalar GFAkid_s2d1=_b[S]

	* Modified Reweighting (S2=1, C=0) 
	qui reg grade_at_age S i.kid_age  [pw=osrweight] if status~=1 & kby<=1997  & kid_age<=15 & kid_age>=6 & ed_flag==1
	return scalar GFAkid_s2d0=_b[S]



	**********************ED by GENDER ***************************************************************************


	*********** Enrollment of own boys**********************************************************************

	* Modified Reweighting (S2=1, C=1) 
	reg enrolled S i.kid_age  [pw=srweight] if status~=1 & kby<=1997  & kid_age<=15 & kid_age>=6 & kid_boy==1
	return scalar ENboy_s2d1=_b[S]


	* Modified Reweighting (S2=1, C=0) 
	reg enrolled S i.kid_age  [pw=osrweight] if status~=1 & kby<=1997  & kid_age<=15 & kid_age>=6 & kid_boy==1 
	return scalar ENboy_s2d0=_b[S]

	*********** Enrollment of own girls**********************************************************************

	* Modified Reweighting (S2=1, C=1) 
	reg enrolled S i.kid_age  [pw=srweight] if status~=1 & kby<=1997  & kid_age<=15 & kid_age>=6 & kid_boy==0
	return scalar ENgirl_s2d1=_b[S]


	* Modified Reweighting (S2=1, C=0) 
	reg enrolled S i.kid_age [pw=osrweight] if status~=1 & kby<=1997  & kid_age<=15 & kid_age>=6 & kid_boy==0 
	return scalar ENgirl_s2d0=_b[S]



	********** grade for age of own boys **********************************************************************
	
	* Modified Reweighting (S2=1, C=1) 
	qui reg grade_at_age S i.kid_age  [pw=srweight] if status~=1 & kby<=1997 & kid_age<=15 & kid_age>=6 & ed_flag==1 & kid_boy==1
	return scalar GFAboy_s2d1=_b[S]

	* Modified Reweighting (S2=1, C=0) 
	qui reg grade_at_age S i.kid_age  [pw=osrweight] if status~=1 & kby<=1997  & kid_age<=15 & kid_age>=6 & ed_flag==1 & kid_boy==1
	return scalar GFAboy_s2d0=_b[S]

	********** grade for age of own girls **********************************************************************
	
	* Modified Reweighting (S2=1, C=1) 
	qui reg grade_at_age S i.kid_age  [pw=srweight] if status~=1 & kby<=1997 & kid_age<=15 & kid_age>=6 & ed_flag==1 & kid_boy==0
	return scalar GFAgirl_s2d1=_b[S]

	* Modified Reweighting (S2=1, C=0) 
	qui reg grade_at_age S i.kid_age  [pw=osrweight] if status~=1 & kby<=1997  & kid_age<=15 & kid_age>=6 & ed_flag==1 & kid_boy==0
	return scalar GFAgirl_s2d0=_b[S]




	restore
end



simulate    kids1=r(kids_s2d1)	kids0=r(kids_s2d0) WHM1=r(WH_M_s2d1) WHM0=r(WH_M_s2d0) HAM1=r(HA_M_s2d1) HAM0=r(HA_M_s2d0) ancM1=r(anemic_M_s2d1) ancM0=r(anemic_M_s2d0)  wp1=r(wp_s2d1)	wp0=r(wp_s2d0)	WHkid1=r(WHkid_s2d1)	WHkid0=r(WHkid_s2d0) HAkid1=r(HAkid_s2d1)	HAkid0=r(HAkid_s2d0) andK1=r(Anemiakid_s2d1) andK0=r(Anemiakid_s2d0) WHboy1=r(WHboy_s2d1)	WHboy0=r(WHboy_s2d0) WHgirl1=r(WHgirl_s2d1)	WHgirl0=r(WHgirl_s2d0)	 HAboy1=r(HAboy_s2d1)	HAboy0=r(HAboy_s2d0) HAgirl1=r(HAgirl_s2d1)	HAgirl0=r(HAgirl_s2d0) 	  ENkid1=r(ENkid_s2d1) ENkid0=r(ENkid_s2d0) 	GFA1=r(GFAkid_s2d1) GFA0=r(GFAkid_s2d0)	 ENboy1=r(ENboy_s2d1)	ENboy0=r(ENboy_s2d0) ENgirl1=r(ENgirl_s2d1)	ENgirl0=r(ENgirl_s2d0) GFAboy1=r(GFAboy_s2d1) GFAboy0=r(GFAboy_s2d0) GFAgirl1=r(GFAgirl_s2d1) GFAgirl0=r(GFAgirl_s2d0)	women_all=r(women_all) women_campaign=r(women_campaign)	, reps(500) seed(8675309): myboot

bstat, stat(non_boot)
